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By a computer simulation approach we study the scattering of p- or s-polarized light from a two- 
dimensional, randomly rough, perfectly conducting surface. The pair of coupled inhomogeneous 
integral equations for two independent tangential components of the magnetic field on the surface 
are converted into matrix equations by the method of moments, which are then solved by the 
biconjugate gradient stabilized method. The solutions are used to calculate the mean differential 
reflection coefficient for given angles of incidence and specified polarizations of the incident and 
scattered fields. The full angular distribution of the intensity of the scattered light is obtained for 
strongly randomly rough surfaces by a rigorous computer simulation approach. 



PACS numbers: 42.25.-p; 41.20.-q 
I. INTRODUCTION 

Theoretical/computational studies of the scattering 
of light from two-dimensional randomly rough perfectly 
conducting surfaces are carried out primarily for two rea- 
sons. These are that a perfectly conducting surface is 
a good approximation to a finitely conducting surface 
in the far infrared region of the optical spectrum, but 
computationally less intensive to study than a finitely 
conducting surface, and that the development of com- 
putational methods for calculations of scattering from 
rough perfectly conducting surfaces can serve as the first 
step in the development of methods that can be used in 
calculations of scattering from rough finitely conducting 
surfaces. 

In the earliest numerical studies of the scattering of 
light from a two-dimensional randomly rough perfectly 
conducting surface [IJ, the pair of coupled inhomogeneous 
integral equations for two independent tangential com- 
ponents of the total magnetic field on the rough surface 
obtained from scattering theory was first converted into 
a pair of coupled inhomogeneous matrix equations by the 
methods of moments . The system of matrix equations 
was then solved by Neumann-Liouville iteration. This is 
a formally exact approach, but one that is computation- 
ally intensive. It is an 0{MN^) approach, where N is 
the number of unknowns to be determined and M is the 
number of iterations 

Subsequent work on this problem has proceeded in two 
directions. One is the exact solution of the integral equa- 
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tions of scattering theory by numerical methods that are 
faster than a straightforward application of the method of 
moments followed by an iterative solution of the resulting 
matrix equation. For example, Wagner et al. j3j have de- 
veloped a fast multipole Fast Fourier Transform method 
to calculate the scattering of an electromagnetic wave 
from a small height two-dimensional randomly rough per- 
fectly conducting surface that is an O(A^lnA^) method. 
For rougher two-dimensional perfectly conducting sur- 
faces they have shown that the multi-level fast multipole 
algorithm, also an 0{Nh\N) method, is more efficient. 

The other direction that has been taken is the approx- 
imate solution of the exact integral equations. In the 
sparse-matrix flat-surface iterative approach of Tsang et 
al. ^E], the matrix elements connecting two close points 
on the surface are treated exactly, while those connect- 
ing two distant points are treated approximately, in an 
iterative solution of the matrix equations obtained by the 
method of moments. This approach has been applied to 
the study of the scattering of electromagnetic waves from 
a two-dimensional randomly rough perfectly conducting 
surface O [7] . It has been elaborated and made faster by 
Johnson and his colleagues, resulting in an 0{N) method 
in some cases, and has been applied to the scattering of 
electromagnetic waves from a two-dimensional randomly 
rough perfectly conducting surface [8J . Soriano and Sail- 
lard |9] have developed a sparse-matrix flat-surface iter- 
ative approach, in which the matrix equations are solved 
by an iterative Krylov method, the biconjugate gradient 
stabilized method [TU] . 

In this paper we return to the approach used in 
where the sparse-matrix flat-surface approximation is not 
used: the matrix elements connecting two points are cal- 
culated accurately for all separations of the two points. 
However, the resulting matrix equations are solved here 
by the biconjugate gradient stabilized method instead of 
by Neumann-Liouville iteration, as in [T]. We show that 
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this approach, together with the increase in computa- 
tional power since [1] was written, provides a simple and 
reliable way of calculating the mean differential reflec- 
tion coefficient for given angles of incidence and specified 
polarizations of the incident and scattered fields, with a 
modest expenditure of CPU time. 

This paper is organized as follows: We start by present- 
ing the scattering geometry considered (Sec. |ll| followed 
by the mathematical formulation of the scattering prob- 
lem (Sec. nil, including the central integral equation on 
which the computer simulations are based. Section IV is 



devoted to the presentation and discussion of the numeri- 
cal results obtained from a rigorous computer simulation 
approach based on an integral equation for the surface 
currents derived in Sec. IIIII A detailed discussion of the 
numerical aspects of such calculations is given in Sec. |V] 
Finally the conclusions that can be drawn from this work 



where 6{z) is the Heaviside unit step function, and 
H(x|a;)i„c is the magnetic component of the incident 
field. 




are presented in Sec. VI 



FIG. 1: (Color online) A sketch of the scattering geometry 
considered in the present work, where the coordinate system 
used and angles of incidence and scattering are defined. 



II. SCATTERING GEOMETRY 

The physical system we consider in this work con- 
sists of vacuum in the region X3 > C(x||), where xy — 
{xi,X2,0), and a perfect conductor in the region X3 < 
C(x||) [Fig. [1]. The surface profile function C(x||) is as- 
sumed to be a single-valued function of xy that is dif- 
ferentiable with respect to and X2, and constitutes a 
stationary, zero-mean, isotropic, Gaussian random pro- 
cess defined by (C(xy)^(X||)) = (5^iy(|xy — Xyj), where 
the angle brackets denote an average over the ensem- 
ble of realizations of the surface profile function, and 
d = (C^(xy))3 is the rms height of the surface. In the 
numerical calculations carried out in the present work we 
will assume a Gaussian form for V[^(|xy — namely 

VFdxy — Xy I) = exp[— (x|| — xj|)^/a^], where a is the trans- 
verse correlation length of the surface roughness. Each 
realization of the surface profile function with these prop- 
erties is generated numerically by a two-dimensional ver- 
sion of the filtering method used in [TT] . 



III. FORMULATION 

A. Integral Equation 

The starting point for our analysis is the Stratton-Chu 
formula [12 for the magnetic field in the vacuum, 

0(a;3-C(xy))H>(x|c^) = H(x|c^)„, 

+ d'x\ [V'.go(x|x')],,^^(,,) X JH(xj||o.), (1) 



The function gQ{x\x') is the scalar free-space Green's 
function and has the representations 



go(x|x') = 



exp [*f|x-x'|] 



(2a) 



(Pq\\ 27rz 



(2^)2 ao(<Zl|) 



exp 



iq|| • (x|| — xji ) 



X exp [iao{qii)\x3 - Xgl] , (2b) 



where 



"o('?||) = ^^"o{q\\) > 0,/TOao((7||) > 0, 



(3) 



and Lo and c are the frequency and speed of light in 
vacuum. In writing Eq. ([T]) we have assumed the 
time dependence exp(— iwi) for the field, but have not 
indicated this explicitly. The (electric) surface cur- 
rent density J//(xy|a;) is defined by J_f/(xy|a;) = [n x 



H>(x|c.)] 



^3=C(X||) 



where n = (-Ci(x||), -C2(x||), 1) is a 



vector that is normal to the surface X3 = C(x||) at each 
point of it, directed into the vacuum, and we have intro- 
duced the notation Cj(x:||) = dC{x\\)/dxj {j — 1,2). 

On evaluating Eq. ([ij at X3 = C(xy) + r] and at 
0:3 = C(xy) — f], where r/ is a positive infinitesimal, adding 
the resulting two equations, and taking the vector cross 
product of the sum with n, we obtain the integral equa- 
tion satisfied by the surface current J//(xy|ci;), 
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Jif(x,i|^) = 23%\^^^\Lu) + Jd^x^nx (|V'.go(x|x')] X JH(xj| 



(4) 



where J^^(x|[|w) = n x H(x|w),„c|^^^^(^||) 
the Cauchy principal value, and we have simplified the 
notation by introducing the definition 



P denotes dition n • Jij(x|||ci;) = 0. Thus only two components of 
Jij(x|||u;) are independent. We choose 3 h{^\\\^)i,2 as 
the independent components, while 



I/(x|x')l - /(x|x') 



2:3=C(X||)' 
^3=C('<:1|) 



(5) 



J_H-(x|||tj)3 = Ci(x||)Jff(x|||tj)i +C2(x||)Jh(x|||cj)2. (6) 



From Eq. Q we find with the aid of Eq. (|6| that 
The system of three equations Q can be reduced to the components J/f (x|| |a;)i_2 satisfy the following pair of 
a system of two equations through the use of the con- equations: 



J«(x||Hi = 2J«(x||Hi--F jd'x\ jfyf (x|||xj|)-,g(°)(x|||xj|)Ci(xj|)-C2(x||).9f (x|||xj|)]jH(xj||c.)i 

.(0)/ 



+5r(x|||x') C2(X||)-C2(X1|) JhKH2 



(7a) 



Jff(x,[|c.)2 = 2J«(x|||c.)2-^P /d^xj, fgr(x|||xj|) Ci(x||) - Ci(xj|) M^[\u;)i 



JO), 



+ 



flf (X|||x') -5^'^^(x|||x')C2(x;|) - Ci(x„)3r(xn|x') Jh{M2 \, (7b) 



.(0) 



where 



5r^(x|||xj,) = |_go(x|x')l-(x,-x;) 



i{uj/c) 

Y _ y'|2 



^'13 



exp[i(a;/c)|x — x'| 



^3=C(X||)' 

4=C(x||) 



(8) 



Equations (|7| are solved by converting them into a pair of coupled matrix equations. This is done by generating a 
realization of the surface profile function on a grid of N"^ points within a square region of the xiX2 plane of edge L, 
where the ratio L/N = Ax is chosen to be Ax = A/7, with A the wavelength of the incident field. The integrals over 
this region in Eqs. are carried out by means of a two-dimensional version of the extended midpoint method [13], 
and the values of J/ffxy |w)i.2 are calculated at the points of this grid. The resulting matrix equations are then solved 
by means of the biconjugate gradient stabilized method |10j . Once J/f (xj[ |cij)i and J//(xj||Lj)2 have been obtained in 
this way, 3 h{^\\\i^)3 is obtained from Eq. ([6]). 

B. Scattered Field 

With the surface current J//(x|||a;) in hand, one is ready to start approaching the calculation of the scattered field. 
To this end, let us start by writing the scattered electric field (in the far zone) in the form 



E(x|tj),c = J j^;^ S{q+,u) cxp[iq+ -x], 



h 

[fp(q+, cj) 7p(q+, cj) + £,.(q+, ^) 7,(q+, w)] exp[iq+ • x], 



(9) 



where = S ■ 'y^ {v = P,s). In writing Eq. Q we have introduced the (unit) polarization vectors ^^(qi, w) 
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for p- and s-polarized scattered light that are mutually 
orthogonal and also orthogonal to the wave- vector q±. 
They can, in accordance with Sipe [H], be defined as 



q± X X3 



qii X X3, 



|q± X X3I 
7s(q±>^) X q± 

g|| X3Tao(g||,^)q|| 



(10a) 



(10b) 



where we have introduced the wave-vector for upward 
(q+) and downward (q-) propagating (plane) waves 



q±(qjht^) = qii ±ao(9ii)x3- 



(lOc) 



From Eqs. (10 1 it is readily shown that the set 
{7p(q±,w), 7jlq±,w), q±(q||,w)} forms a (right- 
handed) orthonormal triad. This implies, for instance, 
suppressing the function arguments for simplicity, that 
7m • 7i^ = ^M-" q± • 7^ = as well as 



7,s = q± X 7p, 

7p = -q± X 7s, 

q± = 7p X 7,. 



(11a) 
(lib) 
(He) 



With the use of one of the Maxwell's equations (Fara- 
day's law), V X E = z(a;/c)H, and Eqs. (Ill, it follows 
from Eq. ^ that the scattered magnetic field can be 
written 



j^:^ [£:p(q+,w)7,(q+,w) - £,(q+, tj) 7p(q+, w)] exp[iq+ -x]. 



(12) 



On the other hand, the scattered magnetic field is also given in terms of the surface current J//(x|| jcj) by the second 
term on the right-hand side of Eq. ([T|, and with the use of Eq. (2b I one is led to (v — p, s) 



2ao(g||) 



Jd^^W 7,.(q+,^) • JH(x|||t^) exp[-iq+ -x]. 



(13) 



The total time-averaged scattered flux is given by the real part of the 3-componcnt of the (complex) Poynting vector 
(S'^ — c/(87r) E X H*) of the scattered field, integrated over the plane 2:3 = 0. From the fields in the form of Eqs. ^ 
and ( 12 1 and the use of Eqs. ( 11 ) we find that it is given by 



P.. 



(14) 



qii is given in terms of the polar and azimuthal scattering angles 9s and (j)s by 



and we recall that q+ = q+(q|| , w), defined in Eq. ( 10c ), depends on the parallel momentum qy . Moreover, the vector 

(15) 



■ sin 6s (cos (jjs , sin 0s , 0) , 



The expression given by Eq. ( 14 1 can then be rewritten as 

c2 1 /wn3 



Ps. 



Sttlu Att^ \ c . 



(^^^ dilsCOs'^Os |£p(q+, w)|^ -I- |£s(q+, w) 



(16) 



where d^ls = sm9sd9sd(f>s is the element of solid angle 
about the scattering direction {9s,(f>s)- 



C. Incident Field 

The incident electric field vector that will be consid- 
ered in this study, has the form of a (Gaussian) beam 



propagating in the direction of 

UJ 

k = — (sin 9o cos 0o, sin 9o sin (f>o, — cos ^o) , (17) 
c 

and is represented by a superposition of incoming plane 
waves 

E(x|w)i„c = J £^'^(q-, w) exp[zq_ • x] 

9||<^ 

X W'(q||,k||), (18a) 
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where M^(qj[,k||) denotes an envelope (or window) func- 
tion, here defined as 



(18b) 



by the plane- wave polarization vectors {e.g. Eqs. (10 1 

and In])), also hold for the polarization amplitudes, £y, 
of the Gaussian beam. 

Moreover, also note that in the limit of a large beam 
width {w oo) Eqs. (20 1 reduce to the plane wave po- 



with w its (and the beam's) half width. Note that in 
the limit of large beam widths {w — > oo), the envelope 
W^(q||,kj[) tends towards (5(q|| — ky) so that, in this limit, 
the incident beam becomes a plane wave. 

A beam as defined by Eqs. (18 1 does not adhere to 



larization vectors given previously in Eqs. ( 10 1 since in 
this limit qn = kii with k\\ = ki. This is another reason 



the usual definition of p- or s-polarized waves since the 
plane of incidence is not well-defined in this case (except 
when w = oo). However, we will still refer to an incident 



beam of the form given by Eqs. ( 18 1 as p-polarized if its 



electric field vector is in the plane "of incidence" defined 
by the vectors k and X3. Therefore, for a p-polarized 
beam, the projection of its amplitude vector £^'^(q_, w) 
onto the xiX2-plane will be parallel to k|| . Moreover, the 

vector amplitude for an s-polarized beam, 5^*''(q_,w), is 
defined as 



for associating the vector amplitudes of Eqs. ( 20 1 with p- 
and s-polarized components, respectively. 

With the polarization vectors available for the incident 
p- and s-polarized components of the incident beam, the 
incident electric field, of given polarization v, can accord- 
ing to Eqs. (18 I and Eqs. (20), be written (assuming unit 
amplitude for simplicity) in the following form 

E^(x|w)j„c = J (fq\\£i'\q-,uj)exp[iq^-x] 



X iy(q||,k|| 



(21) 



In precisely the same way as Eq. (12) was established 



€<i\q-,Lo) = q_ x5«(q_,a;), 



(19) 



similarly to the relation satisfied by the plane- wave po- Eqs. (jUj), that the magnetic component of the incident 



for the scattered field, it follows from Eqs. (21) by us- 
ing Eqs. (19) and relations for 5S>*'' similar to those of 



larization vectors 7^ (c/. Eq. (llal) 



beam then takes the form 



Since in this work we are concerned exclusively with 
isotropic surfaces, we will, with no loss of generality, as- 
sume that the vector ky, if non-zero, is parallel to the xi 
axis, i.e. ky = fcyxi. Under this assumption the ampli- 
tude vector for a p-polarized incident beam, £p\q^,Lu), 
will lie in the a;ia;3-plane, i.e. its second component will 
be zero, which with the condition V • E = (or equiva- 
lently q_ • 5*^*^(q_,a;) = 0) leads us to define 



Hp(x|cj)i„c = J d^gy £i'^(q-,w)exp[iq_ • x] 

xT4/(qy,ky), (22a) 



<}\\<-. 



(20a) 



[9? + ag('7i|)]^ 

The amplitude for the corresponding s-polarized beam 



for a p-polarized beam, and 

Ils{x\uj)inc = - J (fq\\ £^p\q-,Lj)exp[iq^ ■ x] 

X iy(qy,ky), (22b) 



follows from Eq. (19), and, with the use of Eq. (10c 
can be written as 



it 



for an s-polarized beam. 



With the incident field in the form of Eqs. (21) and 



(7192X1 - [qf + al{q\\)]x2 -q2ao{q\\)x3 



(w/c) [qf + a^oin)]' 



([22|, the magnitude of the total time-averaged incident 
fiux is the same for light of both polarizations, and is 
given by 



(20b) 



(23) 



With the beam amplitudes in the form of Eqs. (20) it is 



readily established that similar relations to those satisfied where 



Pine ^ y ^^^911 ao((?y)exp[-'i«^(qy - ky)^] 

= 2ttw'^ (^—^ exp(^-w'^kj^ J dO sin 9 cos'^ 6 Iq (2w'^—k\\ sin 9^ cxp 



2 -2 

-w — sm I 



(24a) 



(24b) 
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r 



and Io{z) is the modified Bessel function of the first kind 
and zero order. In passing, it should be noted that in the 
large beam width limit, for which the beam approaches a 
plane wave, it follows from Eq. (24a I that pinc = 'S'Q;o(fe||) 
where S is the area of the plane — covered by the 
rough surface. 



face, it is the averaged (or mean) of this quantity over 
an ensemble of realizations of the surface that we need 



to calculate. From its definition, we find from Eqs. (14 1 



and (23) that the mean differential reflection coefficient 



for the scattering of incident light of a polarization into 
light of (3 polarization is given by 



D. Mean Differential Reflection Coeflicient 

The differential reflection coefficient is defined as the 
fraction of the total time-averaged flux incident on the 
surface that is scattered into the element of solid angle 
dfls about the scattering direction {9s, (ps)- Since we are 
concerned with scattering from a randomly rough sur- 




If we write the scattering amplitude £^ (q_)_ , uj) as the 
sum of its mean value and the fluctuation about the 
mean. 



each term contributes separately to the mean differential reflection coefficient 



dRfa 



-f-V 

47r2 V c / 



cos 



1 //,in3 

47r2 \ c 



l^/3(q+,'^)r 



1 //,i\3 

47r2 \ c 



l^/3(q+,'^)r 



^/3(q+,'^) 



(26) 



(27a) 



(27b) 



r 



The first term in Eq. ( 27b ) gives the contribution to the 



mean differential refiection coefficient from the light that 
has been scattered coherently. 



dR 



■13a 



dVts 



coh 



— f-V 

47r2 V c / 



cos 



l(g/^(q+,^))r 

Pi7ic 



(28) 



The second term gives the contribution to the mean dif- 
ferential reflection coefficient from the light that has been 
scattered incoherently. 



lm7 



i n coh 



-f-V 

47r2 V c / 



cos^ 9. 



(29) 



The dependencies of the right-hand sides of these ex- 
pressions on the polarization index a is through the de- 
pendence of the amplitudes £/3{q+,uj) on the surface cur- 
rent J^f(x|[|a;) in Eqs. (13 1. This surface current satisfies 
the inhomogeneous integral equations, Eqs. ([7|, in which 
the inhomogeneous terms depend on the incident field, 
and hence on its polarization a = p,s. Thus iJ'^(q+,w) 



depends implicitly on the polarization a of the incident 
field, and so therefore does the differential reflection co- 
efficient. 

The procedure now is to generate a large number Np 
of realizations of the surface profile function C(^||)i ^^id 
for each realization to solve the scattering problem for an 
incident field of p or s polarization. The solution is then 
used to calculate the scattering amplitude 5^(q+, u) and 

|f/3(q+,^)l' 



An arithmetic average of the Np results 



for these quantities yields the q uan ti ties |(£'p(q+, aj))| 
and (|£'s(q+,ti^)P) entering Eqs. (28l-(29l for the mean 
differential reflection coeflicient. 



E. Energy conservation 

To facilitate the discussion of the conservation of en- 
ergy, let us define the following quantity 

U^{9o,^o) = Jdn, (^). (30) 

Recalling the definition of the mean differential reflec- 
tion coefficient, it follows that the physical signiflcance of 
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l4^{6o, 0o) is that it is the fraction of the flux of the inci- 
dent a-polarized light that is scattered into /^-polarized 
light by the rough surface irrespective of scattering di- 
rection. 

For a perfectly conducting surface, all power flux in- 
cident onto the rough surface has to be converted into 
scattered power flux leaving the surface, since there is 
no absorption in the system. Hence, this is nothing but 
energy conservation, and it can be expressed in terms of 

a—p,s 

I3=p,s a=p,s 

= 1, (31) 

where the a-summation over the polarization of the inci- 
dent light is only non-trivial in cases where the incident 
beam does not have a well-defined p- or s-polarization. 
It was pointed out in the previous subsection, that the 
mean differential reflection coefficient can be separated 
into a coherent and an incoherent component. The same 
applies therefore to U{{9o, 0o) and related quantities. 

We note that Eq. pTj ) is rather useful for estimating 
the quality of the simulations, including making sure that 
the discretization interval is flne enough. However, it 
should be stressed that relation (31) is only a necessary 
condition, and its satisfaction does not guarantee that 
the simulations are correct. 



IV. RESULTS AND DISCUSSIONS 

We have carried out calculations of the scattering of 
p- and s-polarized light from a randomly rough perfectly 
conducting surface with an rms height ^ = A and a trans- 
verse correlation length a = 2A, where A is the wave- 
length of the incident field. The polar angles of incidence 
are 6*0 = 0°, 20° and 40°, while the azimuthal angle of 
incidence in all cases is = 0°- The surface is gen- 
erated at a 112 X 112 grid of points covering an area 
= 16A X 16A. The integration mesh size is therefore 
Ax = A/7. The calculations were carried out for an inci- 
dent field in the form of a Gaussian beam [Eqs. ( pO| )] of 
width w = 4A. 

In Fig. |2] we plot the mean differential reflection coef- 
flcients as functions of the polar scattering angle 6^ for 
the in-plane (0^ = 0°) and out-of-plane (0s = ±90°), co- 
(p — > p) and cross-(p — s- s) polarized scattered light due 
to a p-polarized Gaussian beam incident on the surface. 
The results depicted in Figs. [2] were obtained as averages 
over 12, 000 realizations of the surface profile function. 
In obtaining these results we have noted that at least for 
the roughness parameters we have assumed, the contri- 
bution to the mean differential reflection coefficient from 
the light scattered coherently is smaller than the contri- 
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TABLE I: The energy conservation for various polar angles 
of incidence {9o) and incidence polarizations (a) for the sur- 
face parameters given in the text. The surface and scatter- 
ing amplitude were discretized on 112 x 112 and 101 x 101 
grids, respectively. These results were obtained on the basis 
of Eqs. ((sol) and ([31 1. 



bution from the light scattered incoherently by a factor 
of approximately 10"'' (see Table |l] for details). 

There is no single scattering contribution in the cases 
of in-plane cross-polarized [Fig. [2]jb)] and out-of-plane 
co-polarized [Fig. |2jc)] scattering. This we believe is 
the main reason for the reduced amplitude of the mean 
differential reflection coefficients in these cases relative 
to those of Fig. and (d) where single scattering is 
allowed. The peaks at Os = 0° and -20° [20] for in- 
plane co-polarized scattering [Figs. |2ja)] are enhanced 
backscattering peaks [T5J [T71 [T5J HH] . However, the struc- 
tures seen as peaks in the backscattering directions of the 
cross-polarized scattering. Fig. [2jb), are not real peaks, 
as will be seen below from the full angular intensity dis- 
tributions. The results that the mean differential reflec- 
tion coefficients for out-of-plane co- and cross-polarized 
scattering [Figs. |2jc) and (d)] are even functions of 9^ 
are consequences of the scattering geometry, namely that 
(f>o = 0° , (ps = ±90°, and the isotropy of the power spec- 
trum of the surface roughness. 

In Fig. [3] we present corresponding results to those 
of Fig. [3] but now for an s-polarized incident Gaus- 
sian beam. There is no single scattering contribution to 
the in-plane cross-polarized and out-of-plane co-polarized 
scattering, as in the case of p polarization. Also in this 
case the peaks seen in the in-plane co-polarized scatter- 
ing [Fig. (Sja)] are enhanced backscattering peaks, while 
the structures seen in the in-plane cross-polarized scat- 
tering [Fig. |3jb)] in the backscattering direction are not 
real peaks. 

The full angular distribution of the intensity of the 
scattered light is presented as color contour plots in 
Figs. [4}]6| which correspond to the polar angles of in- 
cidence ^0 = 0°, 20°, and 40°, respectively, and for sev- 
eral combinations of the polarizations of the incident and 
scattered light [21]. To the best of our knowledge, this 
is the first time that the full angular distributions of the 
light scattered from a strongly rough surface have been 
obtained by a rigorous computer simulation approach. 
It is observed from Figs. [4[|6|that the angular distribu- 
tions, for given polarizations of the incident and scattered 
light, are far from trivial, and show strong and complex 
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FIG. 2: (Color online) The mean differential reflection coefficients, {dKfiald^s) {a ^ f3), as functions of the polar scattering 
angle Os for the in-plane {4>s = 4>o or (j)s = 4>o + 180°) and out-of-plane {4>s = <^o ± 90°), co- [p — > p) and cross- (p — > s) polarized 
scattering of a p-polarized incident beam (a = p) of width w = AX {6o — 0° and 9o — 20°; 0o = 0°) scattered from a Gaussian 
randomly rough perfectly conducting surface. The Gaussian correlated surface had a correlation length a = 2X and an rms 
height S — X. To facilitate comparison between the various configurations presented in this figure, notice that we have used 
similar scales for all ordinate axes. Moreover, to simplify the presentation of the figures, a convention was adopted where 
negative (positive) values of 6a correspond to <j)s = (jio + 180° {(pa ~ 4'o)- 



angular dependencies. With the full angular dependence 
of the scattered light available, the energy conservation 
of the simulations performed can be obtained by com- 
paring the power incident on the surface to that being 
scattered from it [see Eq. ([3l|)]. For normal incidence, 
we obtained Up — 0.9976 and Us — 0.9970 for p- and 
s-polarized incident light, respectively. For the other an- 
gles of incidence considered, 9q = 20° and 40°, energy 
conservation was satisfied within 0.5% or better (see Ta- 
ble |l] for details) . Even if energy conservation is only a 
necessary requirement, such results, however, still testify 
to the accuracy of the simulations and the approaches 
used to obtain them. 

It is interesting to note that for the roughness param- 
eters considered, the power in a normally incident beam 
is divided essentially equally between p and s polarized 
scattered light (independent of the polarization of the in- 



cident light). This effect we attribute to multiple scatter- 
ing. For the other angles of incidence, it is observed from 
Table |l] that the fraction of incident power being scat- 
tered into the same polarization as that of the incident 
beam (co-polarized scattering), but still independent of 
scattering direction, increases with the polar angle of in- 
cidence. 

We will now discuss Figs. |4j[6]in more detail: We start 
by considering the case of normal incidence; Oq = 0° 
and 00 = 0° [Figs.|4]. Recall that with the assumptions 
and conventions used in this work, the electric field of 
an incident p-polarized Gaussian beam is in the plane of 
incidence. In Fig. |4|a) we present a contour plot of the 
mean differential reflection coefficient for the scattering 
of p-polarized light into either p- or s-polarized scattered 
light, {i.e. the polarization state of the scatted light is 
not being recorded). The angle-dependent scattering, in 
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FIG. 3: Same as Fig. [2j but for an s-polarized incident beam. 



this case, is for the most part rather isotropic, except 
for a sHght anisotropy seen as an elongated (along the 
52-direction) structure around the normal scattering di- 
rection. This structure is caused by the wider intensity 
distribution in the direction perpendicular to the inci- 
dent electric field as compared to the intensity distribu- 
tion along it. The central peak present in Fig. Qa) is 
the enhanced backscattering peak, and is not related to 
specular scattering which for these roughness parame- 
ters can be neglected (see Table |l] for details). A simi- 
lar behavior is seen for the scattering of (normally) inci- 
dent s-polarized light into either p- or s-polarized light 
[Fig. |4|b)]. Here an apparent enhanced backscattering 
peak is also observed. In the case of s-polarization, one 
sees though that the central anisotropic portion of the 
scattering has a different orientation compared to that 
in the case of p-polarization. It remains true, however, 
that there is a stronger scattering perpendicular to the 
(average) direction of the incident electric field indepen- 
dent of the polarization of the incident light. 

Based on these findings, one may be misled into be- 
lieving that the scattering for normal incidence into the 



two possible (linear) polarizations, p or s, is also more-or- 
less isotropic, except maybe for some minor polarization 
dependence for the smaller scattering angles Og- How- 
ever, this is rather far from being true. In Figs.Qc) and 
(d) we present the scattering into p-polarized scattered 
light from, respectively, a (normally) incident p and s po- 
larized Gaussian beam. Similarly, depicted in Figs. Qe) 
and (f ) are the scattering into s-polarized scattered waves 
for an incident p- or s-polarized Gaussian beam. We 
note that taking the sum of the distributions shown in 
e.g. Figs, and (e) produces the angular distribu- 
tion shown in Fig. |4]^a). From Figs. |4][c)-(f) it follows 
that the intensity distributions for scattering from one 
polarization into another, or into the same one, show a 
dipole-like angular dependence. 

For co-polarized scattering, i.e. the polarization of the 
incident light and the (recorded) polarization of the scat- 
tered light are the same, the "forward direction" of the 
dipole-hke pattern is oriented along qi [Figs, ^c) and 
(f)], while for cross-polarization, it is oriented along the 
q2 direction |22j . For normal incidence, the k- vector 
used to define the incident Gaussian beam, does not 
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FIG. 4: (Color online) The complete angular distributions of the mean differential reflection coefficient, {dRpa/ d^s) , for the 
scattering of an a-polarized Gaussian beam incident on the surface at polar angle 9o = 0° and azimuthal angle (/)o = 0°. The 
perfectly conducing rough surface was characterized by a Gaussian height distribution of rms-value 5 = A and a Gaussian 
correlation function of transverse correlation length a — 2A. The incident beam was p-polarized in Figs.|4|a), (c) and (e) [left 
column]; and s-polarized in Figs. [4|b), (d) and (f) [right column]. Moreover, in the top two figures [Figs. [4|a) and (b)] the 
polarization of the scattered light was not recorded; in Figs. [4|c) and (d) [central row] only p-polarized scattered light was 
recorded; while the bottom two figures correspond to recording only s-polarized scattered light [Figs.|4|e) and (f)]. The rough 
surface, covering an area 16A x 16A, was discretized at a grid of 112 x 112 points corresponding to a distcretization interval A/7 
for both directions. The presented figures were obtained by averaging the mean differential reflection coefficient over 12,000 
surface realizations. 
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FIG. 5: (Color online) Same as Figs.[4j but for a polar angle of incidence 9o — 20°. 



(together with X3) define a plane of incidence. How- 
ever, we have used the convention in the simulations, 
that the plane of incidence is defined as the plane having 
(pQ = — sin ^oQli + cos 00^2 as its normal vector which is 
well-defined for all polar angles of incidence (also = 0° ) 
and coincides with the usual definition when 9q ^ Q. 
Since 0o = 0° was assumed for all the simulation re- 
sults presented, it follows (with this convention) that the 



plane of incidence is the gi53-plane. With this definition 
for the plane of incidence, we may rephrase the above 
observation: For co- and cross-polarized scattering the 
dipolc-like pattern is oriented along and perpendicular to 
the plane of incidence, respectively. Later we will see that 
this statement also holds true for non-normal incidence. 

It is noted that we have checked and found that the 
scattering of a normally incident unpolarized beam by 
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FIG. 6: (Color online) Same as Figs.[4j but for a polar angle of incidence 9o — 40°. 



the rough surface, produces, when both its p- and s- 
polarized components are recorded, a fully rotationally 
symmetric intensity distribution (equal to the sum of 
the distributions in Figs. |4|a) and (b)). If only p- or 
s-polarized scattered light is recorded, one will still, with 
the same type of unpolarized illumination, obtain rota- 
tionally symmetric intensity distributions (equal to the 
sum of the distributions from Figs, and (d), in the 



case of p-polarization, and the sum of Figs, ge) and (f) 
for s-polarization) . 

We now turn our attention to the scattering for non- 
normal incidence. In Figs. [5] we present the results for 
the angular distribution of the mean differential reflection 
coefficient for either a p- or s-polarized Gaussian beam 
incident onto the surface at a polar angle 9q = 20° and 
scattered into various polarization states. 
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From Figs, [sj^a) and (b), for which the polarization 
of the scattered light is not recorded, one observes that 
there are pronounced enhanced backscattering peaks lo- 
cated around the backscattering direction (at 9s = 20° 
and (f>s — 180°). It is also observed that the p-polarized 
incident beam tends to scatter more light into the for- 
ward plane {qi > 0) than does an s-polarized incident 
beam. 

The first thing to notice from Figs. |5]^c)-(f), where the 
polarization of the scattered light is recorded, is that the 
co-polarized scattering shows up as an elongated struc- 
ture with the long axis of the pattern directed along the 
plane of incidence, while the cross-polarized scattering 
has the long axis of the scattering pattern perpendicular 
to this plane. This observation is in agreement with what 
was already observed above for normal incidence. How- 
ever, for non-normal incidence, the patterns do show less 
symmetry, as expected, and an even richer and more com- 
plicated angular structure. In principle, the enhanced 
backscattering peak phenomenon should exist in both co- 
and cross-polarized scattering (TTJ [TSl [Hj . However, for 
the roughness parameters assumed in this work, one ob- 
serves instead of a well-pronounced peak in the backscat- 
tering direction, a ridge of constant enhanced intensity in 
parts of the backscattering plane {qi < 0) forming (what 
seems to be) a half circle of constant polar scattering an- 
gle ~ Oq = 20° with e [90°, 270°] [Figs. [Sjd) and 
(e)]. In exactly the backscattering direction, dg = Sq and 
(j)s — 180°, there seems to be little, if any, "extra" en- 
hancement in the cross-polarized scattering as compared 
to the intensities at other values of 0s in the interval 
[90°, 270°]. The enhancement ridge seen is Figs, jsjd) 
and (e) we speculate is caused by a constructive interfer- 
ence effect similar in nature to the underlying enhanced 
backscattering. 

In passing, we note that having available only the in- 
plane and out-of-plane results for the same angle of inci- 
dence, the local enhancements observed in e.g. Figs.[2]^b) 
and[3|b) for 9o = 20°, could easily have been mistaken 
for well-localized features in the backscattering direc- 
tion, similar to what one has for co-polarized scattering 
[Figs.[5](c) and (f)]. In this respect, the angular intensity 
distributions of the kind presented in Figs. |^]6]can pro- 
vide important contributions to a better understanding 
of the multiple scattering phenomena. 

Figures |6] present contour plots of the angular distri- 
butions of the mean differential reflection coefficient for 
a polar angle of incidence 9q — 40°. Since these results 
rather closely resembles those of Figs. [5) we will not dis- 
cuss them further. However, we note that the structures 
due to coherent interference seen in the cross-polarized 
components for 0q — 20°, are much harder to identify in 
the results for Oq = 40°. This is believed to be caused 
by the relatively large angle of incidence, for which it is 
known that coherent effects become weaker |16j . 



V. NUMERICAL ASPECTS 
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TABLE II: The CPU time spent on various stages of the cal- 
culations for one realization of the surface profile function and 
one angle of incidence. All CPU times are measured in sec- 
onds, and the numbers have been rounded to the closest half 
second, and they refer to a machine running an Intel Core2 
CPU (Q9550) operating at 2.83 GHz and running the Linux 
operating system. The surface was discretized on a A'^ x A'^ 
grid of points. The reported CPU times are: the total CPU 
time spent for simulating one surface realization for one an- 
gle of incidence including reading and writing of data (ttot); 
the setup of the system matrix of the linear system Ax — b 
determining the surface currents (tA); the time to solve this 
system by a the iterative BiCGStab method or the direct LU 
decomposition method {tAx=h)', and finally the time to cal- 
culate the refiection amplitudes, £{q+,Lj) for both scattered 
polarizations on a grid of 101 x 101 points (te). The number of 
unknowns to be solved for is A/" = 2A'^^ , where the memory (in 
Gigabytes (Gb)) required to hold the complex system matrix 
A, using single precision, is denoted by Ma oc A/"^ — 4A'^. 

The rigorous computer simulation approach presented 
in this work is rather computationally demanding. 
Therefore, it is important to be able to perform such sim- 
ulations in an efficient manner. One of the most challeng- 
ing aspects of implementing a surface integral method for 
a two-dimensional rough surface, is the memory require- 
ment. By discretizing the relevant integral equations, in 
this case Eq. Q, they are converted into a linear sys- 
tem ^x = b, where A denotes a dense complex system 
matrix; b is the right-hand-side given in terms of the in- 
cident field; and the unknown vector to be solved for, x, 
consists (in our case) of the independent components of 
the surface current J// (xy [a;). If the randomly rough sur- 
face X3 = C(x||) is discretized into NxN points, then the 
number of unknowns would he Af = 2N'^, since for a per- 
fectly conducting rough surface we have two unknowns 
per surface point (the two independent components of 
Jh)- Hence, the amount of memory needed to hold the 
(full) system matrix of the scattering from a perfectly 
conducting surface is A4a = iN^m, where m is the 
size of a single scalar complex variable, which on most 
systems for single and double precision, respectively, is 
ms = 8 bytes and mu — 16 bytes. 

For each surface realization, there are essentially three 
time-consuming steps in this kind of simulation. They 
are: (z) to set up the system matrix; (ii) to solve the lin- 
ear system for the unknown surface currents; and (Hi) to 
calculate the reflection amplitudes. Of the three, it is pri- 
marily the flrst two that are critical and, if not handled 
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properly, particularly the second. For instance, the total 
CPU time taken to complete the calculation using single 
precision and an iterative solver for one angle of incidence 
and one surface realization with N = 112, including read- 
ing input and writing output data, is Uot = 76.0 s on an 
Intel Corc2 CPU (Q9550) operating at 2.83 GHz and run- 
ning the Linux operating system. On the other hand, for 
the same simulation the three steps mentioned above take 
tA — 36.0 s to set up the system matrix, tAx=b = 31.1s 
to solve the linear system by the use of the iterative 
BiCGStab method, and tg = 8.9 s to calculate the re- 
flection amplitudes on a 101 x 101 grid, in total 76.0 s. 
Hence, the additional steps of the calculation, like gen- 
erating the surface, reading and writing data to file etc., 
contribute only insignificantly to the overall CPU time 
{t ~ 0.05s). The computation times for other surface 
discretizations are summarized in Table |lll The reason 
that it takes a relatively long time (compared to ttot) to 
set up the matrix elements is the cost of calculating the 
exponential function contained in the Green's function. 

However, the most critical point to address when try- 
ing to reduce the overall CPU time, is the method used to 
solve the linear system. In this work, an iterative solver 
known as the stabilized bi-conjugated gradient method 
(BiCGStab) [10] has been used, and found to perform 
well and to produce reliable results for our application. 
The iteration process of the BiCGStab solver (using a 
Jacobi preconditioner) was terminated when the relative 
error was 10~^ (or less), which for normal incidence and 
with = 112 required typically a little more then 20 
iterations when starting from an initial guess Xg„ess = 
(of course, other surface parameters and initial guesses 
may require more or fewer iterations in order to reach 
the desired accuracy). Using a direct solver, like the LU- 
decomposition, would have taken significantly longer (see 
Table [ll| . For instance, the time taken to solve the linear 
system for N = 112 by a direct LU solver is 114 times 
longer than that taken by the BiCGStab solver ( Table [ll|. 
Moreover, this difference is expected to increase with in- 
creasing N due to the different scaling with the number 
of unknowns (as also shown by the times presented in 
Table |ll|. It should be noted that a direct solver, like 
the LU-decomposition, opens the possibility for carry- 
ing out calculations for several angles of incidence (the 
right-hand sides of the system) simultaneously with lit- 
tle addition to the overall computation time. This is 
not the case for the BiCGStab-method, where the solu- 
tion time for several angles of incidence scales linearly 
with the number of angles of incidence. There are, how- 
ever, other iterative methods that can solve a linear sys- 
tem with several right-hand-sides with only moderate in- 
crease in computational times. One such method is the 
(restarted) Generalized Minimal Residual Method (GM- 
RES) method mj. Compared to the BiCGStab used 
here, the GMRES is typically more memory demanding 
and, therefore, this possibility has not been explored in 



this work. 

For the sake of comparison, we have repeated the cal- 
culations reported by Tran and Maradudin in Ref. [1] 
using the same numerical parameters (the surface rough- 
ness parameters were already the same) . For the calcula- 
tions carried out in Ref. 1 solving the integral equations 
on a grid of 64 x 64 surface points, each iteration (of 
which there were six) required 365 CPU seconds (on a 
Cray XMP/EA-116 machine), and to calculate the scat- 
tered fields, in-plane or out-of-plane, required 360 CPU 
seconds for each realization of the surface profile func- 
tion, for a total of 2550 CPU seconds for each realization 
of the surface profile function. A similar calculation re- 
quired only 7.6 CPU seconds per surface reahzation, a 
dramatic improvement in speed [23 . This dramatic re- 
duction occurred for two reasons: First, we have the over- 
all improvement in general computer hardware. Second, 
we hold the whole system matrix in memory due to suf- 
ficient memory, while the approach used in Ref. [Ij was 
to regenerate the matrix elements as they were needed. 
This time cost of the latter is not insignificant, as we 
can see from Table [H] and both factors contribute to the 
overall speedup. 



VI. CONCLUSIONS 

In conclusion, we have shown that the use of the 
method of moments and the biconjugate gradient sta- 
bilized method provides a formally exact solution to the 
problem of the scattering of an electromagnetic field from 
a two-dimensional, randomly rough, perfectly conduct- 
ing surface, with a modest expenditure of computational 
time. 

Moreover, the full angular distribution of the intensity 
of the scattered hght, both co- and cross-polarized, was 
obtained by a formally rigorous approach for a strongly 
rough surface. Such distributions can display rather com- 
plex angular patterns that are rooted in the multiple scat- 
tering processes taking place when light interacts with a 
strongly rough surface. 

Due to the full angular intensity distribution being ac- 
cessible, the conservation of energy was checked explicitly 
for all the calculations reported and found to be satisfied 
with an error smaller than 0.5%, or better, something 
that testifies to the accuracy of the approach and a sat- 
isfactory discretization. 
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